home *** CD-ROM | disk | FTP | other *** search
/ Mac Mania 4 / MacMania 4.toast / / Demo's / Igor Demo Pro / 1 PutContentsIn Igor Pro Folder / WaveMetrics Procedures / Analysis / DSP (Fourier etc) / DFTMagPhase < prev    next >
Text File  |  1994-06-03  |  7KB  |  250 lines
  1. #include <BringDestToFront>
  2.  
  3. | Given an input wave (which doesn't require a length that is a power of two),
  4. | create a new wave with the suffix "_dMag" that contains the normalized frequency response
  5. | and optionally create a wave with the suffix "_dPhase".
  6. | The frequency range, and number of output points are specified.
  7. | Options include windowing (Hann) vs no windowing, linear vs dB, phase vs no phase.
  8. | If phase then option of radians or degrees,wrapped or unwrapped.
  9. | You may want to modify the code at the end of this macro.
  10. | It sets the display style and that is a matter of taste.
  11. |
  12. Macro DFTMagPhase(w,window,fMin,fMax,numout,linlog,phase,phasetype)
  13.     string w
  14.     Prompt w,"Input data:",popup WaveList("*",";","")
  15.     variable window=1
  16.     Prompt window,"Windowing:",popup "None;Hann"
  17.     variable fMin=0
  18.     Prompt fMin,"Minimum frequency"
  19.     variable fMax=0
  20.     Prompt fMax,"Maximum frequency, or 0 for Nyquist"
  21.     variable numout=65    | small value recommended.  Remember, we're using the DFT, not the FFT
  22.     Prompt numout,"Number of output points:"
  23.     variable linlog= 2
  24.     Prompt linlog,"Magnitude mode:",popup "Linear;dB"
  25.     Variable phase= 0.5
  26.     Prompt phase,"Phase:",popup "No phase;Phase in radians;Phase in degrees"
  27.     Variable phasetype=1
  28.     Prompt phasetype,"Unwrap phase?",popup,"No;Yes"
  29. ;
  30.     PauseUpdate; silent 1
  31.     
  32.     string destw=w+"_dMag",phasew= w+"_dPhase"
  33.     Variable n= numpnts($w)
  34.     
  35.     if( fMax == 0)
  36.         fMax= 0.5/(pnt2x($w,1)-pnt2x($w,0))
  37.     endif
  38.     if( fMin >= fMax )
  39.         Abort "Min frequency must be less than Max frequency"
  40.     endif
  41.     if( fMax > 1/(pnt2x($w,1)-pnt2x($w,0)) )
  42.         Abort "Maximum frequency too big"
  43.     endif
  44.     if( fMin < 0 )
  45.         Abort "Minimum frequency too small"
  46.     endif
  47.  
  48.     Duplicate/O $w $phasew,$destw
  49.     if(window==2)
  50.         Hanning $phasew; $phasew *= 2            | *= 2 assumes continuous rather than pulsed data
  51.     endif
  52.     Redimension/C $phasew                            |DFT requires complex data
  53.     Redimension/C/N=(numout) $destw
  54.     K0=DFTVarRes($phasew,$destw,fMin,fMax,1)        | Do Forward transform between fMin and fMax
  55.     Setscale/I x,fMin,fMax, $destw
  56.     $destw= r2polar($destw)
  57.     if( phase!=1 )
  58.         Duplicate/O $destw $phasew
  59.         Redimension/R $phasew
  60.         $phasew= imag($destw)
  61.         if( phasetype==2 )
  62.             if( fMin==0)
  63.                 $phasew[0]= $phasew[1]            | try to avoid glitch at dc
  64.             endif
  65.             UnWrap 2*Pi,$phasew
  66.             if( fMin==0)
  67.                 $phasew[0]= 0
  68.             endif
  69.         endif
  70.         if(phase==3)
  71.             $phasew *= 180/Pi
  72.             SetScale y,0,0,"deg",$phasew
  73.         else
  74.             SetScale y,0,0,"rad",$phasew
  75.         endif
  76.     else
  77.         KillWaves $phasew
  78.     endif
  79.     Redimension/R $destw
  80.     if( linlog==2 )
  81.         WaveStats/Q $destw
  82.         $destw= 20*log($destw/V_max)
  83.         SetScale y,0,0,"dB",$destw
  84.     else
  85.         $destw /= n/2
  86.         SetScale y,0,0,"V",$destw
  87.     endif
  88.     BringDestFront(destw)
  89.     if( phase!=1 )
  90.         CheckDisplayed  $phasew
  91.         if( !V_Flag )
  92.             Append/R $phasew
  93.         endif
  94.     endif
  95.     if( numpnts($destw) <= 129 )
  96.         Modify mode($destw)=4,marker($destw)=19,msize($destw)=1
  97.     else
  98.         Modify mode($destw)=0
  99.     endif
  100. end
  101.  
  102.  
  103. | Given an input wave ( doesn't have to be a power of two), compute the real,complex values
  104. | and the magnitude and phase values at the given frequency.
  105. | phase options are radians or degrees
  106. | The results are printed, and are stored in global variables V_real, V_imag, V_mag, and V_phase.
  107. |
  108. | REQUIRES: DFT1()
  109. |
  110. Macro DFTAtOneFrequency(w,freq,phase,printit)
  111.     string w
  112.     Prompt w,"Input data:",popup WaveList("*",";","")
  113.     Variable freq
  114.     Prompt freq,"DFT at frequency (Hz):"
  115.     Variable phase= 1
  116.     Prompt phase,"Phase:",popup "Phase in radians;Phase in degrees"
  117.     Variable printit= 1
  118.     Prompt printit,"Results:",popup "Printed to History;Not Printed"
  119. ;
  120.     PauseUpdate; silent 1
  121.     
  122.     Variable n= numpnts($w)
  123.     String units
  124.     String tmpw=w+"Tmp"
  125.     Duplicate/O $w,$tmpw                            | this could be avoided if we knew whether w was complex
  126.     Redimension/C $tmpw                                | DFT requires complex data
  127.     Variable/C dftCmplx = DFT1($tmpw,1,freq)/(n/2)    | Forward Discrete Fourier Transform
  128.     KillWaves $tmpw
  129.  
  130.     Variable/G V_real= real(dftCmplx)
  131.     Variable/G V_imag= imag(dftCmplx)
  132.     dftCmplx= r2polar(dftCmplx)
  133.     Variable/G V_mag= real(dftCmplx)
  134.     Variable/G V_phase= imag(dftCmplx)
  135.     if(phase==2)
  136.         V_phase *= 180/Pi    
  137.         units= "degrees"
  138.     else
  139.         units= "radians"
  140.     endif
  141.     if(printit==1)
  142.         Print "frequency = ",freq," Hz"
  143.         Print "(magnitude,phase) = (",V_mag,",",V_phase,units,")"
  144.         Print "(real,imaginary)= (",V_real,",",V_imag,")"
  145.     endif
  146. end
  147.  
  148. | Given a complex input wave and a same size output wave, calculate the
  149. | Discrete Fourier Transform using the direct (slow) method.
  150. | parameter dir should be +1 for forward transform or -1 for reverse.
  151. | input wave is not changed
  152. |
  153. | The x scaling of the output wave should be changed by the caller
  154. | with a statement such as
  155. | Setscale/p x,0,1/numpnts(inw)/(pnt2x(inw,1)-pnt2x(inw,0)), outw  
  156. |
  157. | Requires: DFT1()
  158. |
  159. Function DFT(inw,outw,dir)
  160.     wave/c inw,outw
  161.     variable dir
  162.     
  163.     variable n=numpnts(inw)
  164.     variable nn
  165.     variable dx= pnt2x(inw,1)-pnt2x(inw,0)
  166.     
  167.     nn=0
  168.     do
  169.         outw[nn]= DFT1(inw,dir,nn/(n*dx))
  170.         nn+=1
  171.     while(nn<n)
  172.     
  173.     return 0
  174. end
  175.  
  176.  
  177. | Given a complex input wave and a complex output wave (probably of different size),
  178. | calculate the Discrete Fourier Transform using the direct (slow) method for the frequencies
  179. | in the range of fMin to fMax, inclusive.  outw[0] will contain the value for fMin, and outw[n-1]
  180. | will contain the value for fMax.
  181. | parameter dir should be +1 for forward transform or -1 for reverse.
  182. | input wave is not changed.
  183. |
  184. | The x scaling of the output wave should be changed by the caller
  185. | with a statement such as
  186. | Setscale/I x, fMin,fMax, outw
  187. | The values of the output wave have the normal multiplication factor of inPts=numpnts(inw);
  188. |    if the input wave is a record of repetitive data, divide outwave by inPts.  
  189. |
  190. | Requires: DFT1()
  191. |
  192. Function DFTVarRes(inw,outw,fMin,fMax,dir)        | Do transform between fMin and fMax
  193.     wave/c inw,outw
  194.     variable fMin,fMax,dir
  195.     
  196.     variable n=numpnts(outw)
  197.     variable freq
  198.     variable nn
  199.     
  200.     nn=0
  201.     do
  202.         freq= fMin+nn/(n-1)*(fMax-fMin)
  203.         outw[nn]= DFT1(inw,dir,freq)
  204.         nn+=1
  205.     while(nn<n)
  206.         
  207.     return 0
  208. end
  209.  
  210.  
  211. | Given a complex input wave and a same size output wave, calculate the
  212. | Descrete Fourier Transform using the direct (slow) method.
  213. | parameter dir should be +1 for forward transform or -1 for reverse.
  214. | input wave is not changed  
  215. | parameter freq is the frequency at which  the forward DFT will be computed,
  216. | it is the time at which the inverse DFT is computed.
  217. | Returns complex value of Discrete Fourier Transform at given frequency (forward)
  218. | or time (backward)
  219. |
  220. Function/C DFT1(inw,dir,freq)
  221.     wave/c inw
  222.     variable dir,freq
  223.     
  224.     variable n=numpnts(inw)
  225.     variable/c w,wrot
  226.     variable kk
  227.     variable dx=  pnt2x(inw,1)-pnt2x(inw,0)
  228.     variable pointFreq= freq*dx     |pointFreq varies from 0 to 1.  0 for Constant Term, 1 for max frequency term (it plays role of nn/n in DFT Function)
  229.     
  230.     Variable/C out=0
  231.     
  232.     kk=0
  233.     w= cmplx(1,0)
  234.     wrot= cmplx(cos(2*pi*pointFreq),sin(2*pi*pointFreq))    | was     wrot= cmplx(cos(2*pi*nn/n),sin(2*pi*nn/n)) in DFT()
  235.     if(dir== -1)
  236.         wrot= conj(wrot)
  237.     endif
  238.     do
  239.         out += inw[kk]*w
  240.         w *= wrot
  241.         kk+=1
  242.     while(kk<n)
  243.     
  244.     if(dir== -1)
  245.         out /= n
  246.     endif
  247.     
  248.     return out
  249. end
  250.